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ABSTRACT 

Context. Adequate modelling of the multiphase interstellar medium requires optically thin radiative cooling, comprising an inherent 
thermal instability. The size of the occurring condensation and evaporation interfaces is determined by the so-called Field-length, 
which gives the dimension at which the instability is significantly damped by thermal conduction. 

Aims. Our aim is to study the relevance of conduction scale effects in the numerical modelling of a bistable medium and check the 
applicability of conventional and alternative adaptive mesh techniques. 

Methods. The low physical value of the thermal conduction within the ISM defines a multiscale problem, hence promoting the use of 
adaptive meshes. We here introduce a new refinement strategy that applies the Field condition by Koyama & Inutsuka as a refinement 
criterion. The described method is very similar to the Jeans criterion for gravitational instability by Truelove and efficiently allows to 
trace the unstable gas situated at the thermal interfaces. 

Results. We present test computations that demonstrate the greater accuracy of the newly proposed refinement criterion in comparison 
to refinement based on the local density gradient. Apart from its usefulness as a refinement trigger, we do not find evidence in favour 
of the Field criterion as a prerequisite for numerical stability. 
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1. Introduction 

The structure of the turbulent interstellar medium (ISM) of 
the Milky Way and other disk galaxies is determined by 
various complex physical processes. In particular, the for- 
I mation of the neutral atomic phase of the ISM is be- 
lieved to be regulated by condensation under the action of 
thermally bistable radiative cooling. The related thermal in- 
stability (TI) has been extensively studied as a modify- 
ing agent for various driving sources of turbulence in the 
' ISM. Numerical studies include the ma g neto-r otational in- 
' stability (MRI) (Piontek & Ostriker '2005", '2007), explosions 
of SNe (Korpietal. 1999; de Avillez & Breitschwerdt 2005; 
' iJoung & Mac Low 2006; Gressel et al. 2008a b|), and converg- 
, ing flows (I Audit & Hennebellei. 2005; .Vazquez-Semadeni et alj 
120061; iHeitsch et alj l2008l) . Adaptive mesh MHD simulations 
of the last scenario have been performed by i Hennebelle et alj 
(|2008|) . who applied a simple density threshold as a refinement 
criterion. 

The development of TI alone has been explored under vari- 
ous conditions by (Sanchez-Salcedo et al. 2002). Recently, there 
has been a discussion whether TI in conjunction with ther- 
mal conduction can act as an independent source of turbulence 
dBrandenburg et alJl2007l;lKovama & Inutsukall2006h . 

In dynamic simulations of the ISM, the radiative cooling is 
usually treated in the limit of an optically thin plasma, i.e., ra- 
diative transport effects are neglected. The cooling rate is then 
prescribed by A{T) in un its of erg cm-^s~'g~^. H ere we adopt 
the cooling function from lSanchez-Salcedo et all (120 02) that is 
a piecewise power-law fit to results bv lWolfire et alj (I1995D . 



In this paper, we want to focus on an issue in the numerical 
treat ment of therma lly bistable gas that has firstly been brought 
up by lKoyama & Intitsuka (2004 ): The comprehensive theoreti- 
cal analysis by Fieldl (11965 ) reveals that, in the case of vanish- 
ing thermal conduction, TI has a finite growth rate in the limit 
of high wave numbers. This implies that, in the discrete repre- 
sentation of the fluid, numerical fluctuations at grid scale will 
be prone to artificial amplification through the physical instabil- 
ity unless one does take suitable precautions. The most natural 
way to guarantee a physically meaningful, converged solution is 
to introduce a physically motivated dissipation length that then 
needs to be resolved on the numerical g rid. For the c ase when 
thermal conduction cannot be neglected, lFieldl(ll965h derives a 
stability criterion for the condensation mode which has the gen- 
eral form d{pA)/dT - p/T d{pA)/dp < -k^ k with k the co- 
efficient of thermal conduction in units of erg cm"' K"' s"'. In 
the case of a power-law cooling function K{T) ~ this re- 
duces to (jS- l)p^A/r < -k^ K . The wavelength Af where the 
criterion is marginally fulfilled marks the scale where thermal 
co nduction efficiently d amps out fluctuations generated by the 
TI. lKovama & Inutsukai (.2004 ) have shown that 
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has to be resolved with at least three grid cells to attain a con- 
verged numerical solution. It occurs naturally to use this crite- 
rion as a condition for mesh refinement. 

2. Numerical Methods 
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For our simulations, we make use of the newly developed 
version 3 of the NIRVANA code which is a general purpose 
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MHD fluid tool designed to solve multiscale, self-gravitating 
magnetohydrodynamics problems. The NIRVANA code com- 
prises (i) a fully conser vative, divergence-free Godunov-type 
central scheme, Ziegle^ (|2004|) . (ii) block structured adaptive 
mesh refinement (AMR), and (iii) a multigrid-like adaptive mesh 
Poisson solver with elliptic matching, Ziegler (2005). 

The central aim of adaptive mesh simulations is to resolve re- 
gions of interest at e nhanced a ccuracy. In AMR co des like Enzo 
(lO'Shea et allllool. FL ASH (iFrvxell et al.ll2000^ ■ NIRVANA, 
or RAMSES dTevssied l2002h . this can be done by a variety 
of different criteria: baryon or dark matter overdensity (cos- 
mological structure formation), local Jeans-length (protostellar 
core collapse), or entro py gradients (accretion shocks). Recently, 
llapichino et al ] (I2008h have introduced a vorticity-based crite- 
rion to properly resolve a turbulent wake. Alternatively, inde- 
pendent of the underlying physics, the amount of discrepancy 
between the solutions on different AMR levels can be used as an 
indication for refinement CBerger & Colella.l989i) . 

2.1. Classical mesh refinement 

As a standard criterion to trace developing structures, one typ- 
ically monitors local gradients 6u in the fluid variables u. To 
be able to follow moving structures (and to anticipate emerging 
structures), it is desirable to furthermore consider the temporal 
change of these gradients. Due to the additional overhead in stor- 
ing quantities between the timesteps, this is usually discarded. 
Instead, one can make use of the fact that strong gradients are 
bordered by sections of enhanced curvature - reflected in the 
second derivatives 6^u. Accordingly, NIRVANA'S usual mesh re- 
finement is determined by combining the criterion of normalised 
gradients \\6u\\2 : \u\ with one of the ratio of second to first deriva- 
tives ||5^m||2 : ||5m||2 of the conserved variables u, i.e.. 
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where \\6u\\2 and ||(5^m||2 are the first and second numerical deriva- 
tives in the Euclidian norm, and a e [0, 1] controls the mixture 
of the two indicators. A new block is created if the criterion ex- 
ceeds a critical value e„ in a two-cell wide buff'er zone around 
the block to be inserted. Blocks are removed if the criterion falls 
below a factor 0.7 times the critical value. By tuning a, the user 
can specify a bias towards the first or second type of criterion; a 
value of a = corresponds to the standard estimator used in the 
FLASH code jFrvxell et al.ll2000l) . 

For our runs not based on the Field condition, we adopt re- 
finement for the gas density p only and use trigger values €p of 
0.02 and 0.04 while we set a- 0.7. The ratio of grid spacings of 
refinement level I and the base level permits a level-dependent 
refinement via the choice of the exponent ^. Here we fix ^ = 0.6, 
which makes refinements on higher levels successively easier. 
Finally, we set the additional filter value e - 0.01 to prevent the 
refinement of small ripples. 

While this general form of the refinement criterion is very 
flexible, the large number of free parameters renders the opti- 
mal adjustment to a particular problem a tedious task requiring 
a certain amount of experience. Moreover, since AMR always 
implies a trade-off between accuracy and efficiency, the opti- 
mal choice is by no means clear-cut. This is illustrated in Fig.[T] 
where we show the variation of the accuracy and AMR speed- 
ufO] as a function of the two supplementary parameters a and ^ of 
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For a description of the test case, see Sec. l3.1l below. 



Fig. 1. Effect of the adjustable parameters a and ^ on the ac- 
curacy/efficiency of the mesh refinement. Contour lines indicate 
the maximum relative error (in percent, left panel), and the cor- 
responding AMR speed-up (right panel) for the case of ep - 0.02 
and with respect to the fully resolved reference run. 

the classical refinement estimator (|2]l. One can, in principle, use 
these graphs to find the parameter set which yields the highest 
execution speed at a given level of fidelity. However, this proce- 
dure is limited to simple test cases, and its universality cannot 
easily be proofed. 

2.2. Refinement based on the Field condition 

In an alternative approach particularly tailored to trace conden- 
sation interfaces, we apply grid refinement wherever the local 
Field-length Af is not resolved by at least three grid cells. For 
fine tuning, this number could, furthermore, be varied as a lin- 
ear function of the refinement level - although we do not make 
use of this dispensable feature here. The main advantage of the 
proposed new method lies in the fact that it is virtually parame- 
ter free. This implies that the intricate determination of suitable 
values for a and ^ (as described above) can be avoided. 

The evaluation of Af, given by Eq. ([TJ, is straightforward 
since the values of k and T have already been computed for the 
heat flux term. To minimise overheads, they are stored in aux- 
iliary aiTays. It is evident that the definition of the Field-length 
requires the inclusion of a heat conduction term in the energy 
equation. Contrary to this, the FLASH code supports refinement 
based on the local cooling time t - s/p^A, which is indepen- 
dent of a:. If, in our case, k is scaled with p (see below), the two 
methods should perform similarly. 

3. Results 

3. 1 . Gaussian density perturbation 

As a test problem, we revisit the ID case of a den- 
sity perturb ation as c onside red in the model DEN3 from 
Sanchez-Salc edo et al.l (|2002|) . hereafter SSVSG. In this setup, 
the gas with density n - 1 cnT^ initially rests at its equilib- 
rium pressure of p/ks = 2500 Kcm^^* and is perturbed by a 
Gaussian overdensity of amplitude 0.2, FWHM of 3 pc, and 
centred within a domain of 16 pc. Unlike in our simulations, 
SSVSG do not apply thermal conduction. This means that we 
have to choose high enough resolution, such that the conduction 
(which mainly becomes important near the grid scale) does not 
overly affect the features seen in their solution. As a reference, 
we perform a uni-grid run at a linear resolution of 8192 cells, 
corresponding to a grid spacing of ^ 2 x 10"^ pc. The actual 
level of conduction is set such that Af is resolved properly in the 
reference run. This yields a value for k which is about a factor of 
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Fig. 2. Field -length as a function of radius f or the problem DENS 
described inl Sanchez-Salcedo et alj (l2002h for the initial profile 
(dashed lines), and the final profile at f = 4.2 Myr (solid lines). 
The diffusion is set at a level of 1.02 x 10"* erg cm"' K"' s"' and 
scaled as indicated by the labels. Arrows mark a spacing of 3Aa: 
at different refinement levels. 



500 larger than the physical Spitzer value in the limit of vanish- 
ing electron density. The thus even smaller structures expected in 
a physically meaningful scenario illustrate the need for applying 
multiscale techniques. 

Fig. |2] shows a plot of the Field-length for the initial and fi- 
nal density and temperature profiles, where we use three differ- 
ent prescriptions for the conduction: Scaling k with p is typi- 
cally used together with constant kinematic viscosity yielding a 
constant value for the Prandtl-number along with a constant nu- 
merical time step, which is computationally favourable. This is 
also reflected in the fact that the Field-length only shows a weak 
variation with density. For the case of constant k, this variation 
is stronger and accounts for about one order of magnitude in the 
problem under consideration. Finally, for a Spitzer-type conduc- 
tivity with K ~ T^^^, we obtain a variation of over two orders 
of magnitude in Af rendering this case particularly interesting 
for AMR. In the following, we restrict ourselves to the case of a 
constant coefficient /c = 1 .02 x 10"* erg cm"' K"' s"' . 

In general, the results are only marginally affected when con- 
duction is included and agree well with the solution of SSVSG. 
In Fig.[3j the resulting profiles for the run with a 512 cell base- 
grid plus 4 levels of refinement (based on the Field condition) are 
plotted. Before we proceed with an analysis of the AMR perfor- 
mance, we want to compare our reference run 8192-hO (which 
is barely distinguishable from the one shown in Fig.O with the 
profiles from SSVSG: The most noticeable discrepancy is in the 
overall level of the thermal pressure. Due to the assumed peri- 
odicity at the domain boundaries, the mean pressure becomes 
a highly fluctuating quantity (cf. central panel of Fig. [3]). Apart 
from the different offsets, the pressure profiles reasonably match. 
Notably, there is a small pressure dip at the interface of the con- 
densed structure in our runs that is not seen in SSVSG. A similar 
feature is observed in panel (d) of their Fig. 7 (model DEN75), 
although much spikier. This leaves the possibility that the fea- 
ture is not resolved by the plot for model DEN3 in Fig. 2 of 
SSVSG. A comparison run with /<• = in fact shows that, in 
our setup, the dip is broadened by the thermal conduction and 
remains much narrower otherwise. Later investigations on the 
topi c consistently reveal a similar r egion of lower pressure (see 
e.g. IVazquez-Semadeni et al.ll2003h . which can be attributed to 
the fact that the pressure decreases with density in the unstable 
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Fig. 4. Relative errors of the density profiles at time t - 4.2 Myr 
with respect to the fully resolved reference solution. The two 
intervals with the highest deviations are shown for refinement 
based on Af (thick line), Sp - 0.02 (solid fine), and Sp - 0.04 
(dashed fine). 



regime. In this respect, the pressure dip reflects the S-shape of 
the equilibrium cooling curve. 

In Fig. m we plot the relative errors of three AMR runs with 
respect to the fully resolved solution. The two panels correspond 
to the two intervals with the strongest deviations, i.e., an out- 
going wave (upper panel) and the condensation interface of the 
cloud (lower panel). We compare the two runs based on the den- 
sity gradient criterion (with £p - 0.02 and 0.04) with the run 
based on the Field condition. At the time f - 4.2 Myr the num- 
ber of AMR blocks for these runs are 820, 476, and 504, respec- 
tively. Although the number of refined cells is comparable, the 
Field condition yields a relative error that is lower by a factor 
of two at the condensation front. For £p = 0.04 (0.02) the out- 
going wave is resolved at / = 1 (3), while the Field condition, 
expectedly, is insensitive to this feature and does not refine it. 
Still the relative error is lowest in this case, indicating that the 
wave might be overly steepened by the refinement based on the 
density gradients. This finding is rather surprising and deserves 
to be looked at in greater detail. 



4. Discussion 

As is illustrated in Fig.|5] thermal fragmentation produces turbu- 
lent and extremely filamentary structures. Since turbulence can 
only be modelled properly in three dimensions, adaptive mesh 
refinement becomes highly beneficial, if not mandatory. Because 
classical refineme nt strategies reach the ir limits if turbulence is 
involved (see e.g. Ilapichino et al.ll2008l for a vorticity -based re- 
finement criterion), one has to seek for new tracers of structural 
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t = 4.2 




10-1-2-3-4 10 -1 -2 -3-4 10 -1 -2 -3 -4 

log Ixl [pc] log Ixl [pc] log Ixl [pc] 

Fig. 3. Temporal evolution of density (left), thermal press ure (centre), and velocity (rig ht), for our 512+4 run with constant k and 
refinement based on the Field condition (same as Fig. 2 of lSanchez-Salcedo et al ] l200l . 




Fig. 5. The Field-length criterion in action: development of the 
thermal instability from random fluctuations with Spitzer-type 
conductivity. The 2D simulation has a base grid of 256^ with re- 
finement of up to 8 levels, covering scales from 0.003 - 200 pc. 



properties. This is particularly true in the case of a thermally un- 
stable medium. 



(12004^ 



Up to now, the criterion introduced by iKovama & Inutsukal 



is widely disregarded by many members of the nu- 
merical astrophysics com muni t y. Not able exceptions are a 
TI study by Brandenburg et alj ( 20071) and the MRI simula- 
tions of iPiontek & Ostrikerl (l2004l) rThe opposite sta ndpoint, 
repre s ented by a numbe r of authors (see e.g. Gaz ol et alj 
I2005t iJoung & Mac Lowl 12006), is to neglect thermal con- 
duction altogether. The general argumentation is that nu- 
merical diffusion defines a "numerical Field-length" that is 
tho ught to sufficiently damp sm all-scale modes of the instabil- 
ity. [de^wniiZ&SiiticHwerd^ (ilOO^ on the other hand, argue 
that the microscopic heat conduction is suppressed perpendic- 
ular to the magnetic field lines (and thus also isotropically for 
sufficiently tangled fields) and that turbulent transport takes an 
important role. Although this is certainly true, it does not aid the 
discussion of whether explicit conduction should be included. 

In our own simulations, we only observe artificial growth of 
modes near the resolution limit in situations where the Courant 
number is chosen inappropriately high or the gradient based re- 
finement criterion is improperly adjusted. The absence of un- 
physical growth at high wavenumbers might thus, in fact, be 



attributed to the finite level of numerical dissipation, which is 
never really negligible in three-dimensi onal simulations. 

In this respect, the Field criterion of iKovama & Inutsukal in- 
deed seems of academic interest only. Apart from its actual ne- 
cessity, we, however, demonstrate that the very condition can 
successfully be used as a refinement criterion for adaptive mesh 
simulations of the interstellar medium - an approach that is very 
similar to the refin ement procedure b ased on the Jeans-length as 
introduced by Tru elove et al.l (Il997h for self-gravitating clouds. 
In a simple ID test case, the new refinement estimator is found 
to produce more accurate results at comparable numerical cost 
than conventional criteria. 

The overhead of evaluating Ap is low compared to the com- 
putation of the numerical fluxes. Within NIRVANA, the deter- 
mination of the heat flux already requires the evaluation of the 
gas temperature and the (non-uniform) conductivity coefficient. 
These fields can therefore be stored temporarily such that the 
only additional expense comes from evaluating the cooling func- 
tion. 

The definition of Af, nevertheless, implies the inclusion of 
thermal conduction, which is not a standard feature in many 
available codes. An alternative approach (implemented in the 
FLASH code) is to use the cooling time instead, which is in- 
dependent of a:. On the other hand, there is evidence (Piontek, in 
prep.) that heat conduction is in fact physically relevant for the 
formation of molecular clouds, providing further motivation for 
the use of the proposed scheme. 

Compared to the classical mesh refinement based on gra- 
dients, which has to be fine-tuned via multiple parameters ac- 
cording to different situations, the Field criterion is virtually 
parameter-free and gives appropriate results irrespective of the 
particular setup. If, of course, other features like shock fronts 
are supposed to be adequately resolved, one still has to rely on 
a combination with conventional refinement criteria; this does, 
however, not impose any difficulties. 
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